INFO0806 – Statistiques appliquées
FUZELLIER Maxence, BARBET Antoine
24 avril, 2020
Lien GitHub du projet
# Introduction
Dans le cadre du module de statistiques appliquées, nous avions abordé les bases des statistiques et pris en main des outils permettant d’effectuer des operation statistique avancée grâce au langage R, à partir d’un jeu de données. Ainsi, notre choix s’est porté sur des données statistiques relatives à
Ce projet de statistiques appliquées consiste en l’utilisation d’outils R adaptés afin d’effectuer des régression linéaire sur notre jeu de données et produire un rapport.
# Outils et environnement de travail
R et RStudio
R est un langage de programmation dont le but est de pouvoir traiter et organiser des jeux de données afin de pouvoir y appliquer des tests statistiques plus ou moins complexes et se représenter ces données graphiquement à l’aide d’une grande variété de graphiques disponibles. RStudio est une application proposant un environnement de développement et des outils adaptés au langage et à l’environnement de programmation R.
La fonction principale de RStudio consiste à faciliter le développement d’applications en langage R. Pour ce faire, le programme dispose de nombreux outils qui vous permettent notamment de créer des scripts, compiler du code, créer des graphes, ainsi que de travailler avec divers jeux de données.
R markdown
L’extension R markdown permet de générer des documents de manière dynamique en mélangeant texte mis en forme et résultats produits par du code R. Les documents générés peuvent être au format HTML, PDF, Word, et bien d’autres. C’est donc un outil très pratique pour l’exportation, la communication et la diffusion de résultats d’analyse.
Outils complémentaires
- tidyverse : Il s’agit d’une collection d’extensions relatives à la science des données, et permettant la manipulation des tableaux de données, l’import/export de données, la manipulation de variables ou, entre autres, la visualisation de données ;
- data.table : Extension et manipulation avancée des tableaux de données ;
- plotly et highcharter : Comme ggplot2 (compris dans le tidyverse), ces packages permettent la visualisation de données, à la différence près qu’ils rendent ces graphiques interactifs ;
- hrbrthemes : Thèmes pour ggplot2.
Données
Chargement des jeux de données :
data <- read_csv("data/prestige.csv") %>% as_tibble()
data## # A tibble: 102 x 6
## education income women prestige census type
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 13.1 12351 11.2 68.8 1113 prof
## 2 12.3 25879 4.02 69.1 1130 prof
## 3 12.8 9271 15.7 63.4 1171 prof
## 4 11.4 8865 9.11 56.8 1175 prof
## 5 14.6 8403 11.7 73.5 2111 prof
## 6 15.6 11030 5.13 77.6 2113 prof
## 7 15.1 8258 25.6 72.6 2133 prof
## 8 15.4 14163 2.69 78.1 2141 prof
## 9 14.5 11377 1.03 73.1 2143 prof
## 10 14.6 11023 0.94 68.8 2153 prof
## # … with 92 more rowsRégression linéaire
En statistiques, et en apprentissage automatique, un modèle de régression linéaire est un modèle de régression qui cherche à établir une relation linéaire entre une variable, dite expliquée, et une ou plusieurs variables, dites explicatives.
Ici, nous testons donc la relation entre la variable prestige (il s’agit d’un score de prestige relatif à la profession) et la variable éducation (qui reflète le niveau d’étude).
education <- data %>%
ggplot(aes(x = education, y = prestige)) +
geom_point() +
geom_smooth(method = loess, color = "red", fill = "#69b3a2", se = T) +
geom_smooth(method = lm, color = "blue", se = F) +
theme_ipsum()
scatter_edu <- ggplotly(education)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'scatter_eduNous testons aussi la relation entre la variable prestige et la variable income.
income <- data %>%
ggplot(aes(x = income, y = prestige)) +
geom_point() +
geom_smooth(method = loess, color = "red", fill = "#69b3a2", se = T) +
geom_smooth(method = lm, color = "blue", se = F) +
theme_ipsum()
scatter_inc <- ggplotly(income)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'scatter_incprest.lm1 <- lm(prestige~education, data = Prestige)
prest.lm1##
## Call:
## lm(formula = prestige ~ education, data = Prestige)
##
## Coefficients:
## (Intercept) education
## -10.732 5.361bacf <- acf(residuals(prest.lm1), plot = F)
bacfdf <- with(bacf, data.frame(lag, acf))
conf.level <- 0.95
ciline <- qnorm((1 - conf.level)/2)/sqrt(length(residuals(prest.lm1)))
lag_plot <- bacfdf %>%
ggplot(aes(x = lag, y = acf)) +
geom_hline(aes(yintercept = 0)) +
geom_segment(aes(xend = lag, yend = 0)) +
geom_hline(aes(yintercept = ciline), linetype = 3, color = 'darkblue') +
geom_hline(aes(yintercept = -ciline), linetype = 3, color = 'darkblue') +
theme_ipsum()
lag_plot <- ggplotly(lag_plot, tooltip = c("lag", "acf"))
lag_plotdurbinWatsonTest(prest.lm1)## lag Autocorrelation D-W Statistic p-value
## 1 0.2752512 1.43917 0.006
## Alternative hypothesis: rho != 0qq_edu <- prest.lm1 %>%
ggplot(aes(sample = residuals(prest.lm1) / 10))+
stat_qq() + stat_qq_line() +
#theme_ipsum() +
labs(x = "Theoretical Quantiles\nlm(prestige ~ education)",
y = "Standardized residuals")
#qq_edu <- ggplotly(qq_edu)
qq_edu
shapiro.test(residuals(prest.lm1))##
## Shapiro-Wilk normality test
##
## data: residuals(prest.lm1)
## W = 0.98065, p-value = 0.1406prest.lm2 <- lm(prestige~income, data = Prestige)
prest.lm2##
## Call:
## lm(formula = prestige ~ income, data = Prestige)
##
## Coefficients:
## (Intercept) income
## 27.141176 0.002897qq_inc <- prest.lm2 %>%
ggplot(aes(sample = residuals(prest.lm2) / 10))+
stat_qq() + stat_qq_line() +
#theme_ipsum() +
labs(x = "Theoretical Quantiles\nlm(prestige ~ income)",
y = "Standardized residuals")
#qq_inc <- ggplotly(qq_inc)
qq_inc
shapiro.test(residuals(prest.lm2))##
## Shapiro-Wilk normality test
##
## data: residuals(prest.lm2)
## W = 0.97169, p-value = 0.02729plot(prest.lm1, 3)
ncvTest(prest.lm1)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.6327545, Df = 1, p = 0.42635plot(prest.lm1,1)
summary(prest.lm1)##
## Call:
## lm(formula = prestige ~ education, data = Prestige)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.0397 -6.5228 0.6611 6.7430 18.1636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.732 3.677 -2.919 0.00434 **
## education 5.361 0.332 16.148 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared: 0.7228, Adjusted R-squared: 0.72
## F-statistic: 260.8 on 1 and 100 DF, p-value: < 2.2e-16confint(prest.lm1)## 2.5 % 97.5 %
## (Intercept) -18.027220 -3.436744
## education 4.702223 6.019533influenceIndexPlot(prest.lm1)
outlierTest(prest.lm1)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## newsboys -2.98896 0.0035306 0.36012prest.lm1bis <- lm(prestige~education, data=Prestige[-c(53,67),])
compareCoefs(prest.lm1 ,prest.lm1bis) ## Calls:
## 1: lm(formula = prestige ~ education, data = Prestige)
## 2: lm(formula = prestige ~ education, data = Prestige[-c(53, 67), ])
##
## Model 1 Model 2
## (Intercept) -10.73 -11.26
## SE 3.68 3.53
##
## education 5.361 5.417
## SE 0.332 0.318
## my_df <- data.frame(education=c(10.25))
predict(prest.lm1, newdata=my_df)## 1
## 44.21701predict(prest.lm1, newdata=my_df, interval="prediction")## fit lwr upr
## 1 44.21701 26.06518 62.36885predict(prest.lm1, newdata=my_df, interval="confidence")## fit lwr upr
## 1 44.21701 42.40008 46.03395my_df <- data.frame(education=c(10.25, 11.25, 12.25))
predict(prest.lm1, newdata=my_df,interval="confidence")## fit lwr upr
## 1 44.21701 42.40008 46.03395
## 2 49.57789 47.75810 51.39768
## 3 54.93877 52.89190 56.98564my_pres <- Prestige
my_pres$res <-residuals(prest.lm1)
head(my_pres)## education income women prestige census type res
## gov.administrators 13.11 12351 11.16 68.8 1113 prof 9.250875
## general.managers 12.26 25879 4.02 69.1 1130 prof 14.107621
## accountants 12.77 9271 15.70 63.4 1171 prof 5.673573
## purchasing.officers 11.42 8865 9.11 56.8 1175 prof 6.310758
## chemists 14.62 8403 11.68 73.5 2111 prof 5.855950
## physicists 15.64 11030 5.13 77.6 2113 prof 4.487854my_pres$fitted <-fitted(prest.lm1)
head(my_pres)## education income women prestige census type res
## gov.administrators 13.11 12351 11.16 68.8 1113 prof 9.250875
## general.managers 12.26 25879 4.02 69.1 1130 prof 14.107621
## accountants 12.77 9271 15.70 63.4 1171 prof 5.673573
## purchasing.officers 11.42 8865 9.11 56.8 1175 prof 6.310758
## chemists 14.62 8403 11.68 73.5 2111 prof 5.855950
## physicists 15.64 11030 5.13 77.6 2113 prof 4.487854
## fitted
## gov.administrators 59.54913
## general.managers 54.99238
## accountants 57.72643
## purchasing.officers 50.48924
## chemists 67.64405
## physicists 73.11215head(predict(prest.lm1))## gov.administrators general.managers accountants purchasing.officers
## 59.54913 54.99238 57.72643 50.48924
## chemists physicists
## 67.64405 73.11215head(fitted(prest.lm1))## gov.administrators general.managers accountants purchasing.officers
## 59.54913 54.99238 57.72643 50.48924
## chemists physicists
## 67.64405 73.11215vcov(prest.lm1)## (Intercept) education
## (Intercept) 13.520978 -1.1835055
## education -1.183505 0.1102162ggplot(Prestige, aes(y=prestige, x=education))+
geom_point()+
geom_smooth(colour="red", method="lm", fill="red") +
ylab("Prestige")+
xlab("education") +
theme_classic()+
annotate("text", x = 9, y = 80, label = "prestige = -10.73 + 5.36 * education\n (pval<0.001)")## `geom_smooth()` using formula 'y ~ x'
int_pred <- predict(prest.lm1, interval="prediction")## Warning in predict.lm(prest.lm1, interval = "prediction"): predictions on current data refer to _future_ responsesmy_prest2 <-cbind(Prestige, int_pred)
head(my_prest2)## education income women prestige census type fit
## gov.administrators 13.11 12351 11.16 68.8 1113 prof 59.54913
## general.managers 12.26 25879 4.02 69.1 1130 prof 54.99238
## accountants 12.77 9271 15.70 63.4 1171 prof 57.72643
## purchasing.officers 11.42 8865 9.11 56.8 1175 prof 50.48924
## chemists 14.62 8403 11.68 73.5 2111 prof 67.64405
## physicists 15.64 11030 5.13 77.6 2113 prof 73.11215
## lwr upr
## gov.administrators 41.33302 77.76523
## general.managers 36.81573 73.16903
## accountants 39.52816 75.92469
## purchasing.officers 32.33470 68.64379
## chemists 49.31584 85.97226
## physicists 54.67820 91.54609ggplot(my_prest2, aes(y=prestige, x=education))+
geom_point()+
geom_smooth(colour="red", method="lm", fill="red") +
geom_line(aes(y=lwr), color = "red", linetype = "dashed")+
geom_line(aes(y=upr), color = "red", linetype = "dashed")+
ylab("Prestige")+
xlab("education") +
theme_classic()+
annotate("text", x = 9, y = 80, label = "prestige = -10.73 + 5.36 * education\n (pval<0.001)")## `geom_smooth()` using formula 'y ~ x'